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Heterogeneities in Aging Models of Granular Compaction 
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Kinedcally constrained models (KCM) are systems with trivial thermodynamics but often complex dynamical 
behavior due to constraints on the accessible paths followed by the system. Exploring these properties, the Kob- 
Andersen (KA) model was introduced to study the slow dynamics of glass forming liquids and later extended 
to granular materials. In this last context, we present new results on the heterogeneous character of both in and 
out of equilibrium dynamics, further stretching the granular-glass analogy. 
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I. INTRODUCTION 

In recent years, the analogy between structural glasses and 
dense granular systems has become deeper and extensively 
explored filE ipl i S i i E 

[17, 18, 1 91 Ha mi m. As temperature is lowered, or den- 
sity increased, respectively, both systems undergo a glass, or 
jamming, transition, where the relaxation times dramatically 
increase. Yet, in spite of all evidence of phenomenological 
similarity, the two systems are fundamentally different, for ex- 
ample, in the length scale of its components (molecules versus 
macroscopic particles) and the role of thermal energy (none, 
in the case of granulars). Because the thermal energy is too 
small to induce movement in these macroscopic particles, en- 
ergy should be externally supplied, for example, by vibrating, 
tapping, shearing, rotating, etc, the system, in order for a gran- 
ular system explore its configurational space. 

Recently, the role of dynamical heterogeneities in the com- 
plex dynamics of glassy systems has been addressed f23, 
|24, 25, 26]. The glass transition seems to be purely dy- 
namical, with no increasing static correlation length as the 
system approaches the transition, differently from the usual 
critical slowing down. Thus, the increase in relaxation 
times seems to be related to a diverging dynamic correlation 
length 1 27. .2&. .29„ „30> .3 L 32. .33,1 , associated to the increas- 
ing number of particles whose displacements become dynam- 
ically correlated, the heterogeneities, that develop during the 
evolution of the system. In Ref. 1 34] we showed, at the level 
of a simple, kinetically constrained model issll . that dynami- 
cal heterogeneities seem to play the same role in granular sys- 
tems as they do in structural glasses, with an increasing length 
scale as the system approaches the jamming transition. This 
has also been observed in other recent works I10i.i36i.l37ll3^ . 
However, the precise role played by these structures and the 
associated lengths, on the dynamics of granular and colloidal 
systems, is yet to be understood. 

The jamming transition has also been studied in con- 
nection with the concept of dynamically available volume 
(DAV) L39.. .40.1 . Apart from having enough nearby empty 
space I4lll . a particle is mobile if, for the particular type of 



model considered here, the displacement is allowed by the ki- 
netic constraints. Empty sites that are able to receive a neigh- 
boring mobile particle are called holes. Holes can be classi- 
fied as either connected or not, the former being those that, 
by allowing a particle to jump into it may eventually facili- 
tate the movements of all particles in the system. On the other 
hand, nonconnected holes leave a backbone of blocked parti- 
cles. Close to dynamical aiTest, the density of connected holes 
decreases with the density of particles and is related with the 
inverse of the bootstrap length Ii42..t43il . the average distance 
between two connected holes, in bootstrap percolation. In 
turn, this length is associated with the transport properties of 
the system L44„ .45, .46J . 

Here we make an initial attempt of studying how these holes 
behave when the system is externally driven and falls out of 
equilibrium, in particular, their role during the compaction 
regime of granular systems, what is connected with the voids 
distribution measured experimentally |47]. Although we do 
not distinguish at this stage connected from disconnected 
holes, this would be important in order to fully understand 
the microscopic compaction mechanism. 



II. KOB-ANDERSEN MODEL 
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The Kob-Andersen model 114 ST is one of the simplest mod- 
els describing the complex dynamics of glassy and granular 
systems. It consists of a lattice gas of particles, each site 
being either empty or occupied by one particle, with no static 
interactions between them, i.e., y{ —0. In addition, a kinetic 
constraint |35] should be satisfied in order to allow the dis- 
placement of a particle to an empty neighboring site: there 
should be fewer than m occupied nearest neighbors before 
and after the move. This kinetic rule is time-reversible and 
detailed balance is satisfied. If the constraint is obeyed, the 
particle is said to be mobile and the companion vacant site is 
said to be a hole. At high densities, the dynamics slows down 
because the reduced free-volume makes it harder for a parti- 
cle to satisfy the dynamic constraints. These constraints were 
introduced to mimic the cage effect where, due to geometri- 
cal effects, the displacements of the particles are hindered by 
their neighbor particles. 

Although on hypercubic lattices this model presents a jam- 
ming transition only at full occupancy |44], with a super Ar- 
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rhenius behavior of the relaxation time, for finite lattices the 
bootstrap length may become larger than the system size L. 
Thus, for a finite sample, the system may be frozen due to the 
lack of connected holes, implying the existence of a dynam- 
ical critical density Pc(i-) < 1, slowly increasing with L. At 
this point, it is observed that the diffusivity falls to zero as a 
power law |48], £>(p) ^ (Pc — p)*^, with p^ depending, in addi- 
tion, both on the lattice geometry and on the kinetic constraint 
m f34Ll43.l49ll . Assuming such a power law form for the diffu- 
sivity, finite systems, both with and without gravity, are well 
described, qualitatively and quantitatively, by a nonlinear dif- 
fusion equation |34, 50, 51]. 

Since much of dynamical properties of both structural 
glasses and dense granular materials are dictated by steric con- 
straints, we have generalized 1 52] the Kob-Andersen model by 
including a gravitational field. The Hamiltonian now has a one 
body term, = jY.iZini, where n,- = 0, 1 is the occupation 
variable of the i-th site whose height is z,-, y = mg/k^T is the 
inverse gravitational length and g is the constant gravitational 
field acting in the downward direction. We follow a contin- 
uous vibration dynamics, assuming that the random diffusive 
motion of particles, produced by the mechanical vibrations 
of the box, can be modeled as a thermal bath of temperature 
T. The particles satisfying the kinetic constraints may always 
move downwards while upward movements are accepted with 
a probability x ~ exp(— y), related to the vibration amplitude. 
Particles are confined in a closed box of bcc structure, with 
periodic boundary conditions in the horizontal direction. We 
set the constraint threshold at m = 5. As the Markov process 
generated by the kinetic rules is irreducible on the full config- 
uration space is^l . the static properties of the model are those 
of a lattice-gas of non-interacting particles in a gravitational 
field, and these can be easily computed. 

III. HETEROGENEITIES 

Several measures for quantifying spatial heterogeneities 
have been introduced for kinetic models |35]. In particular, 
this issue was recently investigated in the KA model 145. jBT 
without gravity using fourth-order correlation functions i42l 
I43I1 . In Ref . i34il these were extended to include the non-zero 
gravity case. In Fig. [2 we plot the dynamical nonlinear re- 
sponse 

X4{z,t)=N(^{q\zj))-{q{z,t))^) (1) 

where is the number of sites in the computation, q{z,t) — 
C(z,f)/C(z,0) and 

CfcO = ^L«<(fK(0)-p(z,f)pfcO). (2) 

For all vibrations considered, the system is in the fluid phase 
and is able to achieve the asymptotic state very fast. Notice 
that, differently from previous works, here we measure these 
quantities separately around each layer of the system (to im- 
prove the averages, / runs over all sites in the z, z — 1 and z+l 
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FIG. 1: Dynamical response, eq.Q as a function of time (in MCS) 
for different vibrations: x=0.92 (filled symbols) and 0.94 (empty 
symbols). Different symbols stand for different heights: z = 5, 10 and 
15 (square, circle , tri angle, respectively). The line is the asymptotic 
X4 = 1 behavior Notice in the figure the presence of two very 
close curves: they correspond to different vibrations and heights, but 
their density is the same to within numerical accuracy. From I34ll . 

layers). In other words, we are probing horizontal, micro- 
scopic heterogeneities, not the macroscopic ones due to the in- 
trinsicly inhomogeneous vertical profile. Analogously to what 
happens in the KA model without gravity and in other glassy 
systems, the peak is shifted to higher times and gets larger as 
the density increases (the lower is z, the greater is the density), 
which is an indication of cooperative dynamics, as larger clus- 
ters have more difficulty to respond to a perturbation. As these 
measures are done for vibrations above the apparent jamming 
threshold, the behavior of the system should not be affected by 
the finite size shortcomings discussed above and one is able 
to obtain the true, infinite size behavior Indeed, the growth of 
the peak is compatible with the known relaxation time L44.1 . 

T expexp (^Y3^) ^'^^ 

as can be seen in Fig. |2] Analogous results can be obtained 
also from the peak of the Kovacs hump |54] and persistence 
times distribution. Interestingly, X4 only depends on z through 
its local density, t4{z^t') — X4(p(z),f)- Fig.[l]shows that two 
curves corresponding to different heights and vibrations (z and 
X), but having almost the same density (within numerical pre- 
cision), coincide. 

A more direct measure of inhomogeneity comes from the 
mobile particles or, equivalently, the holes at different heights, 
without distinguishing, at this stage, between connected and 
non-connected holes. In the inset of Fig.|3 we show the den- 
sity and holes profiles, p(z) and v(z), respectively, still in the 
fluid phase, as a function of z for two different vibrations x. 
The density decreases with height and the stronger the vibra- 
tion, the flatter is the density profile and the broader is v(z), 
that is, the more vertically homogeneous the system is. In- 
deed, for X — \ (no preferential direction) both profiles are 
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FIG. 2: Position of the peak in X4- Tiie growth is compatible with 
the known relaxation time i44ll . Eq.|3|(ciashed line). 



flat. Notice also that as x decreases, the holes concentrate in 
the region that will form the interface once the system goes 
out of equilibrium. The interesting result that can be seen in 
Fig.|3lappears when we eliminate z and plot v(z) as a function 
of p(z): the data for different values of x perfectly collapse 
onto a single master curve. Moreover, they are on top of the 
points for the KA model without gravity. In other words: even 
if the system is no longer homogeneous, the local density of 
holes in equilibrium only depends on the local density, and not 
on the whole profile, or on the vibration or the height. This, 
together with the analogous result for X4, if general, is an im- 
portant property: as some of the quantities depend only on the 
local density, not on the whole profile, there is further sup- 
port for even simpler, one dimensional models as well as local 
density approximations |34, 50, 51]. >From the point of view 
of simulation, it is a fast way to obtain, in just one simulation, 
the whole profile for systems without gravity. 

Upon decreasing the vibration, a finite system enters in the 
aging regime and the density profile develops two very dis- 
tinct regions SsdR : an almost flat plateau at Pc for z< zo and a 
density decreasing region (interface) for z > zo- The position 
of this plateau is size dependent and slowly increases with L, 
similar to p^. Analogously, the hole profile also has two cor- 
responding regions (see Fig.l^and the inset): at the interface 
where most holes are localized and in the bulk, where their 
number is much smaller These regions are separated by a 
dip with a few layers width, where there are almost no holes, 
corresponding to the dense layer seen in the density profile 
near the interface, at the top of the granular pile (this region 
also seems to become more localized with time), as seen in 
Fig. 13 At the interface, the profile is strongly peaked, and 
moves to the left, accompanying the interface, as the system 
ages (compactifies). On the other hand, in the bulk, the profile 
is linear, with a tangent that slowly decreases in time (although 
also compatible with the inverse logarithm law, a good fit is 
obtained with f^"-^^). This evaporation of holes is the direct 
mechanism of the compactification process. Also, because of 
the almost hole-free layers between the bulk and the interface. 




0.2 0.4 0.6 0.8 1 
P 

FIG. 3: Density of holes as a function of density for x = 0.92 (trian- 
gle) and 0.94 (circle), along with the results for the KA model with- 
out gravity (line) in a bcc lattice of height 4L, with L = 20. Inset: 
holes profile (symbols) along with the corresponding density profiles 
(lines). The higher x corresponds to the flatter density profile and to 
the broader holes profile. 



it is very hard to exchange particles between the bulk and the 
interface. The interesting result here is that this process is in- 
homogeneous: it is faster at the topmost region of the bulk 
(just below the dense layer) and slower at the bottom of the 
system |55|, as seen in Fig.|5] The oddity comes from the 
fact that the compaction is faster where the density is higher, 
the region where one would expect a slower evolution and a 
smaller number of holes. Instead, the evolution is faster and 
the number of holes is higher. Finally, it should be remarked 
that the small positive gradient seen in the density profile is 
consistent with experimental results f?,*?], although its sign is 
model dependent [ 56.1 . 

In Fig. |5]the parametric plot of v versus p in the out-of- 
equilibrium, aging regime is shown for several vibrations x 
along with the equilibrium curve. Notice that despite the fact 
of not being stationary, the points, for several times (but the 
same x), roughly fall on a master curve, while different vi- 
brations no longer collapse onto the same curve. However, 
these x-dependent universal curves no longer correspond to 
the equilibrium one, and as the vibration increases, there is 
a drift toward the equilibrium curve. The two regions seen 
in the v profile, Fig.0] contribute differently for the curve in 
Fig. |6l the bulk forms the high p, small v region near the 
horizontal axis (enlarged in the inset), while the interface gen- 
erates the rest of the curve. In the inset of Fig.|6]we can see 
that the bulk behavior is the opposite of what would be ex- 
pected at high densities: as the density increases the corre- 
sponding holes also increase! This explains why the density 
profile evolves faster near the dense layer where the density is 
higher: the greater the number of holes, the easier it is to com- 
pactify. Moreover, for different times, the curves no longer 
collapse (notice that the inset of Fig.|S]shows data for differ- 
ent times but same vibration x). Thus, differently from the 
holes from the interface, there is no universal curve for the 
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FIG. 4: Hole profile at different times (smaller time at the top). The 
profile is made of two clearly separated parts: ;) the interface (inset) 
contains almost all holes in the system and moves to the left as the 
system becomes compact; /;) the bulk has a much smaller hole frac- 
tion, that increases linearly with height. Separating these regions we 
have a few layers near z = 100 that are almost hole-free. The de- 
clivity slowly decreases with time. Notice also the difference on the 
vertical scales. 
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FIG. 5: Bulk density profile at different times (smaller time at the 
bottom). The region that compactifies the most is the denser one, 
close to the top, due to a larger number of holes. 



holes in the bulk. Interestingly, besides having the opposite 
p dependence, for larger times the curves move away from 
the equilibrium one (that is one order of magnitude above). 
Even the interface behavior is puzzling: layers whose density 
is small would be expected to be in local equilibrium. How- 
ever, the hole distribution only coincides with the equihbrium 
one for p — » 0. 




FIG. 6: Density of holes as a function of density for several values 
of X in the aging regime (0.2, circle, and 0.4, triangle), along with 
the equilibrium results for the KA model (line). Differently from the 
fluid phase, here there is no collapse of the curves onto the equilib- 
rium curve. Notice that there are some points with p > Pc that come 
from the oscillatory/dense layer regions. It is interesting to notice 
that in the case without gravity, the number of holes in the high-p re- 
gion is below the equilibrium one, while here it is above. Inset: Holes 
coming from the bulk for x = 0.2 and three different times (from left 
to right: t = 1072, 10974 and 99999) as a function of density. Dif- 
ferently from what would be expected, this an increasing function 
of density. In these plots we do not take into account the oscillating 
layers near the bottom and near the dense layer, only the bulk layers 
(6 < z < 95). 



IV. CONCLUSIONS 

In summary, besides briefly reviewing the first application 
of the dynamical susceptibility X4 in the context of granular 
compaction |34|, we also extended the notion of dynamically 
available volumes for this externally driven, out of equilib- 
rium situation. This is an additional similarity between struc- 
tural glasses and granular systems, while shedding some light 
on the microscopic mechanisms responsible for the slow dy- 
namics close to these transitions. Interestingly, the stationary 
behavior of the Kob-Andersen model with gravity is charac- 
terized by the local density, in spite of the macroscopically in- 
homogeneous density profile and the vibration imposed: lay- 
ers with the same density presents the the same time depen- 
dence of X4{z,t) as well as the same density of holes. In the 
aging, compaction regime, this is no longer the case, although 
there is still a certain degree of universality in the behavior 
of the total number of holes from the interface, reflected on 
the collapse, at different times (but different densities), of all 
points onto a roughly universal curve that depends, in its turn, 
on the vibration (temperature). However, how connected and 
non-connected holes contribute to these profiles are yet to be 
investigated. 

In addition, there are still many issues that deserve a closer 
inspection. Important information can also be obtained from 
the different sectors of the X4(0 function |28]. However, the 
range of time/density considered here is still small to resolve 
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these different sectors. The distinction between connected and 
disconnected holes should also be made, as the former should 
be much more important in the compaction mechanism, along 
with the study of their spatial distribution |57]. Correlations 
may be also defined considering only the holes, with the asso- 
ciated dynamical susceptibility. 
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